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Ultrahigh energy cosmic rays from 
collisional annihilation revisited 

R. Dick, K.M. Hopp, and K.E. Wunderle 



Abstract: We re-examine collisional annihilation of superheavy dark matter particles in dark 
matter density spikes in the galactic halo as a possible source of ultrahigh energy cosmic 
rays. We estimate the possible flux in a way that does not depend on detailed assumptions 
about the density profiles of dark matter clumps. The result confirms that collisional 
annihilation is compatible with annihilation cross sections below the unitarity bounds for 
superheavy dark matter if the particles can form dense cores in dark matter substructure, and 
it provides estimates for core sizes and densities. The ensuing clumpy source distribution 
in the galactic halo will be tested within a few years of operation of the Pierre Auger 
observatory. 
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1. Introduction 

The source of the observed ultrahigh energy cosmic rays (UHECRs) with energies beyond the 
Greisen-Zatsepin-Kuzmin cutoff, E > Egzk — 4 x 10^^ eV, remains a puzzle for several reasons: 

1. Scattering off cosmic microwave background photons limits the penetration depths of charged par- 
ticles at these energies to distances < 100 Mpc [1-3]; 

2. the distribution of arrival directions of UHECRs does not seem to favor any known astrophysical 
sources within the GZK cutoff length; 

3. it seems extremely difficult to devise sufficiently efficient astrophysical acceleration mechanisms 
which could accelerate particles to energies E > Egkz', 

see e.g. [4-6] for recent reviews. For the second point, we note that Dolag et al. [7] find that typical 
deflection angles due to magnetic fields should remain below the angular resolution of the cosmic ray 
observatories for UHECRs originating within 100 Mpc, but the pattern of arrival directions does not 
seem to point to any known astrophysical accelerators in that range (but see [8] for a different opin- 
ion on the possible size of deflection angles in bottom-up models). It has been pointed out that there 
may be a correlation with BL Lac objects at cosmological distances [9]. However, generation in BL 
Lacs would require both an extremely powerful acceleration mechanism and conversion into a neutral 
component to avoid the GZK cutoff. This difficulty is avoided in "local" bottom-up scenarios like e.g. 
UHECR generation in local gamma-ray bursts within a radius of 90 Mpc [10], but Scully and Stecker 
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had pointed out that gamma-ray bursts will not provide the observed flux [11]. Proposals of non-local 
explanations for absence of a GZK cutoff include a possible violation of Lorentz invariance at high 
energies, see [12] for a recent discussion. A recent examination of angular correlations by Gorbunov 
and Troitsky confirmed the absence of correlations with nearby visible objects [13]. 

The difficulties to explain extremely high energetic cosmic rays had motivated Hill already in 1983 
to propose monopolonium decay as a possible source [14]. More recently the constraints 1-3 had also 
motivated Berezinsky, KachelrieB and Vilenkin [15] as well as Kuzmin and Rubakov [16] to propose 
decay of superheavy dark matter particles as sources of the observed UHECRs, see also [17]. The 
terrestrial flux should then be dominated by decays in our galactic dark matter halo. Decay of super- 
heavy dark matter is an attractive proposal because it explains UHECRs from simple conversion of 
rest mass energy in our cosmic neighborhood, thus evading the GZK bound and the need for extremely 
powerful astrophysical sources. However, it is difficult to identify decay mechanisms of particles of 
mass Mx > 10^^ GeV which are slow enough to ensure a lifetime tx > 10^" years. This had moti- 
vated the authors of Refs. [18-20] to propose coUisional annihilation of free superheavy dark matter 
particles (SHDM particles, wimpzillas) as another mechanism to generate UHECRs in our cosmic 
neighborhood. 

The lifetime problem is avoided in the annihilation scenario because the conversion of rest mass 
energy into UHECRs proceeds through two-particle collisions, thus providing a kinematic suppression 
of the reaction rate fix oc — n^(c7Af )• The three UHECR puzzles are solved as in the decay scenar- 
ios: direct conversion of rest mass energy eliminates the need for an extremely powerful acceleration 
mechanism; the UHECRs originate within our galactic halo; and the distribution of arrival directions 
is determined by dense dark matter subclumps in the galactic halo, rather than by the locations e.g. of 
active galactic nuclei. 

However, the coUisional annihilation scenario suffers from its own specific shortcoming: the minute- 
ness of the annihilation cross section a ~ of supermassive particles on the one hand ensures that 
these particles can survive to the present day (even with about equal amounts of matter and anti-matter), 
but on the other hand requires high-density substructure in the galactic halo to generate the observed 
flux. The results of [18] used two separate models for dark matter subhalo profiles and indicated a 
strong dependence of the flux on the shape of the dark matter clumps. This left open whether col- 
hsional annihilation could be a generic possibility or might be constrained by strong restrictions on 
subhalo profiles. We present a simplified estimate of the flux from coUisional annihilation in SHDM 
clumps to answer this question. The main purpose of our simplified estimate is to clarify the depen- 
dence of the expected flux from SHDM annihilation on key parameters like core densities and core 
radii of dark matter subhalos in the galactic halo. We find that the hypothesis of coUisional SHDM 
annihilation is compatible with the observed flux if the cores of dark matter subclumps of the galactic 
halo can reach densities of order 10~^ qq and extensions of a few AU. This corresponds to a mass 

IO^'Mq and would indicate that about 10% of the mass of the smaUest (and most abundant) sub- 
structures detectable with A/'-body simulations would be bound in dense cores^. We denote these dense 
cores as "wimpzilla stars". 

Besides the Ufetime constraint for relic decay models and the cross section constraint for annihila- 
tion, top-down scenarios in general may be constrained by indications that the highest energy cosmic 
rays could be hadron dominated, see e.g. [23]. On the other hand, careful re-analysis showed that the 
observed highest energy event could have arisen from a primary photon [24]. 

Photon domination would be expected e.g. if superheavy neutrinos annihilate through an interme- 
diate Z or a Higgs boson, because the ensuing tt^ shower would decay into UHE photons over mean 
free paths of order 10 km. The composition of showers in such a model would therefore resemble Z- 



^ Observational support for the existence of CDM subhalos within galactic halos comes from gravitational lensing 
[21,22]. 
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bursts, which assume annihilation of ultrahigh energy neutrinos (instead of supermassive neutrinos) 
through scattering with background neutrinos [25-27]. However, the different origin of Z-bursts im- 
pUes a markedly different anisotropy signature: Z-bursts at reUc neutrinos in our GZK volume would 
produce an isotropic pattern of arrival directions without connection to the galactic halo. 

The first crucial experimental test for our top-down model of SHDM annihilation in WIMPZILLA 
stars will be the observed anisotropy signature after a few years of operation of the Pierre Auger 
observatory [28]: 

If our annihilation model is correct, a pointlike source distribution should emerge with increasing 
source density towards the galactic center [20]. In particular, the majority of observed UHE events 
after about four years should appear in doublets or higher multiplets. 

Alternatively, if the decay model is correct, the anisotropy signature after about five years should trace 
the galactic dark matter halo with uniformly increasing intensity towards the galactic center. 

In Sec. 2 we discuss the generation of weakly coupled particles during inflation and explain why 
stable particles with masses m > 10^^ GeV are still around to create the observed UHE cosmic rays 
through annihilation. In Sees. 3 and 4 we simplify the estimates from [18] for the UHE cosmic ray flux 
which might be expected from SHDM annihilation in dark matter substructures within our galactic 
halo. The method developed in Sees. 3 and 4 yields improved estimates of typical core sizes and 
densities of these substructures. Sec. 3 reviews the basic formalism for the calculation of cosmic ray 
fluxes from dark matter annihilation, and in Sec. 4 we present our simplified estimate of the flux from 
SHDM annihilation. We set ?i = c = 1 except in equations which are used for numerical evaluations. 
The differential flux per energy interval is j{E) = dl/dE, and the flux per unit of soUd angle is 
J{E) = (Pl/dEdn. 

2. Superheavy dark matter from non-ad iabatic expansion 

One of the most interesting sources for superheavy dark matter particles is gravitational production 
during non-adiabatic expansion. This mechanism is usually considered in terms of the Bogolubov 
transformation between in and out vacua [29, 30], and was also extensively discussed as a source of 
superheavy dark matter particles [31-38]. Here we would Uke to point out how this mechanism for 
particle production can be understood in terms of the evolution equation for weakly coupled helicity 
states in an expanding universe. In the present paper we will discuss this only for weakly coupled 
helicity states which are not ultra-relativistic (A ^ m~^), although we are sure that similar conclusions 
can be drawn for the ultra-relativistic case, too. 

The direct couphng between non-adiabacity and changing particle numbers during gravitational 
expansion follows directiy from the Einstein equation and the cosmological principle, since it is well 
known that the evolution equations for the Robertson- Walker metric 

a^ + k K 
and 

„a d'^ + k 

2- + 5— = -Kp 

imply energy balance with only mechanical energy: 

d{ga^) = -pda^, (2) 

and therefore the first law of thermodynamics spUts into Eq. (2) and a separate balance equation for 
entropy and particle densities: 

Td{sa^) = - ^ Hid{via^). (3) 
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At the same time Eq. (1) usually implies a decrease or at best constancy of the comoving energy 
density a^{t)Q{t) during expansion. Consider e.g. the standard Unear dispersion relation for radiation 
or dust 

P=g-l)., (4) 

where the parameter is £ — 4 for radiation and £ — 3 for dust. After inflation, when k can be neglected 
and radiation initially dominates, the scale factor and dominant energy density evolve according to 
a{t) oc i^/^ and g{t) oc t^^, leading to a decrease a^{t)g{t) oc of the comoving energy density. 
However, during inflation we have p :^ —q leading to constant energy density q and an exponential 
increase a^{t)Q{t) oc exp(3i?i) of the comoving energy density. From this argument we expect an 
energy increase at least in all those states which cannot readily redistribute their energy through suffi- 
ciently strong interactions, and inspection of the corresponding evolution equations for weakly coupled 
helicity states confirms this: 

Towards the end of inflation the curvature parameter k is negigible on subhorizon scales, and weakly 
coupled helicity states have to satisfy 

^(k, t) + 3H4>{\si, t) + (m^ + exp(-2ilt)) 0(k, t) = 0. (5) 

Here (j){}s., t) is defined through Fourier expansion with respect to the dimensionless comoving coordi- 
nates X 



(l>{x,t) = 3 / (i'^k(/)(k, t) exp(iaik • x). 

•\/27r J 



27r 

The factor ai = a{ti)exp{—Hti) is related to the scale factor a{ti) at the beginning of inflation. 
Towards the end of inflation the exponential factor in Eq. (5) is exponentially small, exp{—2Ht) « 
exp(— 140), and we find the same approximate evolution equation for all modes which are not ultra- 
relativistic (|k| ^ m) at the end of inflation: 

(f>{\si, t) + 3H^{\si, t) + mV(k, t) ~ 0. (6) 

H is constant during inflation, and therefore the time evolution of the comoving energy density in the 
weakly coupled helicity states is given by 

a'{t)Q^{\i,t) ~ ^a3(i) (,^(k,t)<^(-k,t)+mV(k,t)</'(-k,i)) 

~ A+ exp(t\/9H^ - 4m2) + A_ exp(-t\/9H'^ - Am'^^ + B. 

This implies a growing mode in the comoving energy density of weakly coupled states with m < 
1.5H ~ 10^'* GeV. What is special about the superheavy particles is that their comoving energy density 
is conserved after inflation, because the behavior of massive (m > t~^) weakly coupled helicity states 
in the subsequent radiation and dust dominated backgrounds preserves their energy: 
The asymptotic solution for weakly coupled massive states with m> in such a background yields 

(/){]<., t) (X t~^l^ cos(mi -h ^p) , 

£»0(k,t) oc oc a~^, 

and this implies in particular that the comoving density of massive particles freezes out at the end of 
inflation (i ~ lO'^^ s) if 

m > i"^ ~ 10^^ GeV. 
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From these considerations follows a mass window for inflationary production of superheavy rehc par- 
ticles 

10^2 GeV < m < lO^GeV. 

It is certainly reassuring for the top-down models that a typical time scale for the duration of inflation 
yields the UHECR mass scale for the preserved particles. 

3. Basic formalism 

Following [18], we consider annihilation of two superheavy dark matter particles of mass Mx > 
10^^ GeV primarily into two jets of energy Mx- 

The cosmic ray flux per energy interval originating from annihilation of particles and anti-particles 
of total density nx(r) is 

if the particles are not Majorana particles, or 4 x j{E, Vq) otherwise. dM {E, Mx) is the number of 
particles in the energy interval [E, E + dE] emerging from a jet of energy Mx and two primary jets 
are assumed. Neglecting events with more than two jets, dhf {E, Mx)/dE is related to fragmentation 
functions via 

dE ^^2aAdE~2Mx^ l^^^^^^xj, 

i i 

where x = E/Mx and F'-^^ (x, 4M|^) is the differential number of particles of species i generated in 
the prescribed x-range in an annihilation event with s = 4M|-. 

As shown in [18], UHECRs from annihilation of superheavy dark matter particles provide already 
a good fit to the shape of the observed spectrum if the old MLLA fragmentation function [39] is 
used for an estimate on dM{E, Mx)/dE. More refined approximations have been employed in the 
discussion of the decay scenario [17,40^3], and these also further improve the corresponding fit for 
the annihilation scenario, since for the spectral shape the transition between the decay and annihilation 
scenarios only corresponds to a scaling of jet energies 



dAf{E,Mx/2) 



dAf{E,Mx) 



decay 



dE 

in the fragmentation function 



(8) 

annihilation 



4. A simplified estimate of the flux from dark matter clumps in the galactic 
halo 

Numerical simulations indicate that the dark matter halos of galaxies have a lot of substructure in 
terms of spatial density fluctuations. In [18] it was found that the flux from local overdensities, or dark 
matter clumps, exceeds the flux from the smooth halo component by far. Here we will reconsider these 
calculations with a view on the density nx and size rcore of the central core regions of dark matter 
clumps, nx is of crucial importance for the UHECR flux from annihilation in these clumps. 

We also want to simplify the flux estimate with regard to assumptions on the distribution of sub- 
clumps. In [18] a clump distribution 

-3/2 

(9) 
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was assumed, where Uci {d, Md) is the volume density of clumps with mass Md at distance d from the 
galactic center, do is a scale radius for the subclump distribution in the galaxy, and Mh is the mass of 
the galactic halo [44J. Some of the recent work on structure formation has questioned the assumption 
of increasing clump density towards the galactic center because tidal stripping would more severely 
affect large clumps close to the galactic center. However, the distribution of the 150 known globular 
clusters in the galaxy suggests an increase of visible substruture towards the galactic center [45], and 
it seems reasonable to assume a similar distribution for dark matter substructure. 

In the present paper, we want to avoid any detailed assumptions like (9) about the profile of the 
substructure distribution in the galactic halo, and we also want to avoid any assumptions about the 
density profiles of substructure, except for the reasonable assumption that the substructures will have 
dense cores with an average mass density nxMx- 

The annihilation flux from an approximately compact source of volume V at distance d{'> V^^^) 

is 

VP ^ [ ^3 1 dAf{E,Mx) , ,2/ X V dN{E,Mx) 2 , \ 

^(^''•^) = y/^ 8.|re-rp dE ""^^^ " 8^ dE ^^'^ 

UHECR experiments often plot the approximately flat rescaled flux 

vi:. ^ dJ\f{x,Mx) 2 / V 

^ " S^dUlxc^ d. 

where x = E/Mxc^ is energy in units of primary jet energy. A typical value for fragmentation func- 
tions at a; = 0.1 is dAf{x, Ejet)/dx ~ 30. 

In the following we parametrize the annihilation cross section in terms of the s-wave unitarity 
bound for Mc^ lO^^ GeV and v = lOOkm/s: 

{aAv) = e X ^ = ^ X 4.40 X lO'^^ ^^3/^^ (11) 

where the unitarity bound is ^ < 1. 

For a fiducial average distance we use d = 10 kpc, whereas for the fiducial volume and density of 
subclump cores we should have 

nxMxVcore ^ O.lMd = Q.l.fclMhalo/Nd ~ IQ-H'Ihalo- 

Here fd is the mass fraction in the N^i dark matter subclumps. Numerical A^-body simulations and 
gravitational lensing indicate that fd should be of the order of a few percent. E.g. a core density 

nxMx = VQe = V ^ 1-41 x lO'^kg/m^, 

= ^ ^ 10i4eV/c^ = ^ ^-^^ ^ 
yields typical core volumes and radii 

Vcore ^ T]-^ X2X 10^ Vq, 
fcore ^ V~^^^ X 1-9 X lO" m. 

With these fiducial average values, the flux at E = 10^^ GeV from Nd dense cores of clumps can be 
estimated as 



E-^m NdV^ore ^ I' 2'^^'''^''^ n\{aAv) (12) 

£=1011 GeV S,TTd^MxC^ dx 

E^ dM{x,Mx) 

©2008 NRC Canada 



R. Dick, K.M. Hopp, K.E. Wunderle 7 
^ 5.83 X lO^'^eV^m-^s-i x M x . . x ( x 



0.06 2 X IQI^Mq VlOkpcy IQ-^' 
E'^JiE) ~ 4.64 X lO^^eV^m-^s-^sr-i (13) 

£=1011 GeV 

fci Mhalo f d r?e 

X-—— X X — — X 



0.06 2 X 1O12M0 VlOkpcy 10-3 ' 
This is compatible with the UHECR flux observed by AGASA^ [46], 

E^J{E)\E=iQiiGeY - 4.2 X 102^eV^m-2s-isr-\ 

if the product of the parameters rj and ^ is of order r]^ ~ 10"^. E.g. reproduction of the beautiful 
standard model fragmentation fit in Fig. 7 from Ref . [41] with the present approximations would require 
7?^ ~ 8.6 X 10-4. 



5. Conclusion 

The result (12) indicates that coUisional annihilation of superheavy dark matter particles remains 
a viable scenario for the origin of UHECRs, albeit at the expense of postulating dense cores in dark 
matter substructure. This difficulty compares e.g. to the lifetime problem in decay scenarios, or to the 
problem of defeating synchrotron radiation loss and coUisional energy loss in bottom-up scenarios. 

The unitarity bound (1 1) on the annihilation cross section is both a virtue and a curse of the coUi- 
sional annihilation scenario. It is a virtue because it implies that ultrahigh energy cosmic rays from col- 
Usional annihilation must originate in dense cores of dark matter clumps in the galactic halo, amounting 
to a unique prediction: an unmistakable clumpy pattern in arrival directions should emerge after four to 
five years of operation of the Pierre Auger observatory, appearing as pointlike sources within our halo. 
The connection to the halo should become apparent through an increase of source density towards the 
galactic center, and it will eliminate extragalactic scenarios for UHE origin. In addition, the clustering 
of arrival directions distinguishes it from the SHDM decay scenarios, and the absence of correlations 
of the pointUke sources with supemovae will rule out a bottom-up origin. A further advantage of the 
unitarity bound is that it ensures that stable superheavy dark matter particles will still be around due to 
the extremely slow rates of any reactions involving these particles and whatever affected the low mass 
modes to generate the baryon asymmetry after inflation very likely could not affect the superheavy 
modes. Yet at the same time the unitarity bound seems to require high core densities, and this may 
colhsional annihilation look less appeahng in the current stage, where it is not exphcitly ruled in or 
ruled out. However, the required core densities for dark matter subclumps are not unreasonable and the 
possibility of UHECRs from coUisional SHDM annihilation can not easily be dismissed. 

From our point of view, primary composition is a secondary criterion for separation of differ- 
ent scenarios for UHECR origin. Inferring primary composition from observations requires particle 
physics modelling with extrapolations to very high energies whereas the interpretation of anisotropy 
will involve much less theoretical uncertainty. 

Either way, the anisotropy signature observed by the Pierre Auger observatory will decide the fate 
of the annihilation scenario and we are cautiously optimistic that the coUisional SHDM annihilation 
mechanism in WIMPZILLA stars will withstand the tests of the near future. 



^The HiRes collaboration uses a different binning of the AGASA data, but still finds a similar AGASA flux at and 
above E = 10^^ GeV without a cutoff, see e.g. [47]. However, the HiRes collaboration reports a cutoff in the 
spectrum based on its own data. This issue will be resolved by the Pierre Auger Observatory due to the much 
improved statistics expected from that experiment, and due to the combination of Cherenkov detectors and air 
fluorescence telescopes. 
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